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Driven Diffusive Systems: How Steady States Depend on Dynamics 
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In contrast to equilibrium systems, non-equilibrium steady states depend explicitly on the un- 
derlying dynamics. Using Monte Carlo simulations with Metropolis, Glauber and heat bath rates, 
we illustrate this expectation for an Ising lattice gas, driven far from equilibrium by an 'electric' 
field. While heat bath and Glauber rates generate essentially identical data for structure factors 
and two-point correlations, Metropolis rates give noticeably weaker correlations, as if the 'effective' 
temperature were higher in the latter case. We also measure energy histograms and define a simple 
ratio which is exactly known and closely related to the Boltzmann factor for the equilibrium case. 
For the driven system, the ratio probes a thermodynamic derivative which is found to be dependent 
on dynamics. 
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I. INTRODUCTION 



In statistical physics, Monte Carlo (MC) simulations 
play a major role for the study of phase transitions and 
critical phenomena, as well as ordered and disordered 
phases []|. Leaving out many details of the art of com- 
puting, the broad outline of the simulation process is eas- 
ily summarized. For systems both in and far from ther- 
mal equilibrium, dynamic rates W[a — > a'] are defined 
which specify how a given configuration a evolves into a 
new one, a' , when an update is attempted. The simula- 
tions then generate long sequences of such configurations. 
Once initial transients have decayed, time-independent 
(stationary) observables - which will be our focus in the 
following - can be computed as configurational averages. 
For a system in thermal equilibrium, characterized by 
a Hamiltonian TL, it is well known that any choice of 
W's, as long as they satisfy detailed balance with respect 
to H, will generate configurations distributed accord- 
ing to the same Boltzmann factor, exp(—{3H). In other 
words, time-independent observables, including both uni- 
versal and non-universal properties, are independent of 
the choice of rates, provided detailed balance holds. The 
resulting freedom can be exploited to designparticularly 
efficient codes, such as cluster algorithms 2]. In stark 
contrast, no such 'decoupling' of dynamic and stationary 
characteristics occurs for systems driven out of equilib- 
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rium: even though non-equilibrium steady states (NESS) 
display time- independent observables, these are sensitive 
to modifications of the dynamic rates. This behavior can 
be traced back to the violation of detailed balance which 
is an inherent feature of non-equilibrium systems 0, 0] . 

While the sensitivity of NESS to the choice of the dy- 
namics has been noted before 0, Q , no systematic com- 
putational study has yet been undertaken. In this pa- 
per, we consider two models: the standard Ising lattice 
gas and its non-equilibrium cousin, the driven Ising lat- 
tice gas (or KLS model, after the initials of its inventors 
Both involve particles diffusing on a lattice, sub- 
ject to an excluded volume constraint and an (attractive) 
nearest-neighbor interaction. The total number of par- 
ticles remains conserved. For the Ising lattice gas, the 
rates for particle hops to unoccupied nearest-neighbor 
sites are chosen to satisfy detailed balance, with respect 
to the Ising Hamiltonian. In contrast, the driven version 
involves an additional external force which acts on the 
particles much like an electric field on (positive) charges: 
aligned with a lattice axis (e.g., y), it favors particle hops 
along its direction. In conjunction with periodic bound- 
ary conditions, this bias breaks detailed balance and es- 
tablishes a non- equilibrium steady state. This NESS dif- 
fers drastically from its equilibrium counterpart, exhibit- 
ing generic long-range correlations, a novel universality 
class, and highly anisotropic ordered phases 0. 

To illustrate the importance of the transition rates for 
the NESS, we measure structure factors and two-point 
spatial correlations in the driven case, using Metropolis 
0, Glauber and heat bath j^] rates. To date, simula- 
tions of the driven model have focused almost exclusively 
on Metropolis rates other rates have only been in- 
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voked in some analytic studies [1113. While the first two 
are easily implemented for conserved particle number, an 
appropriate generalization of heat bath rates is designed 
here for the first time. To avoid complications due to 
inhomogeneous ordered phases, we choose temperatures 
above or at criticality. For comparison, we also show 
the same quantities for the (equilibrium) Ising model: as 
expected, they are found to be identical up to statisti- 
cal fluctuations, and almost perfectly isotropic. In stark 
contrast, when the drive is turned on, we find strongly 
anisotropic behavior and markedly different values for 
the three choices of rates. While all driven cases exhibit 
the signatures of long-ranged decay in real space, the cor- 
relations are much weaker for Metropolis rates than for 
either Glauber or heat bath ones. 

As an additional probe into the differences between 
the standard Ising lattice gas and its driven cousin, we 
construct energy histograms (with respect to the Ising 
Hamiltonian) for both. For the equilibrium case, these 
are of course intimately related to the Boltzmann distri- 
bution and contain a wealth of thermodynamic informa- 
tion. For the driven system, they are easily measured 
in a simulation, but their physical interpretation has not 
been established yet. Here, we consider a simple his- 
togram ratio whose equilibrium limit is easily derived, 
and compute its non-equilibrium counterpart. 

This paper is organized as follows. We first introduce 
our models and the three types of dynamics, followed by 
a brief discussion of our key observables. We then present 
our data for two-point correlations, structure factors and 
energy histograms. We conclude with a summary, offer- 
ing a conjecture for the origin of the observed differences 
between the three rate functions. 



II. BACKGROUND 



A. The models 

In this section, we introduce our two prototype models, 
namely, the Ising lattice gas and its driven version. Both 
are defined on an M x L square lattice in two dimensions, 
with fully periodic boundary conditions. Each site i is 
either occupied by a particle or empty, which we denote 
by a spin variable o~\ taking two values: +1 (occupied) 
or —1 (empty). For the equilibrium Ising model, we can 
specify a (global) Hamiltonian: 



H[a] 



-J 



J2 ai °h 

<i,j> 



(1) 



where the sum runs over nearest- neighbor pairs of sites, 
and J > denotes the binding energy. In order to access 
the Ising critical point, we consider only half-filled sys- 
tems: (LAI )~ 1 v\ — 0- When coupled to a heat bath 
at temperature T, the probability, Pq(o~), to find the sys- 
tem in configuration a is controlled by the well-known 



Boltzmann factor: Pq(ct) cx exp(— (3H) with f3 = l/ksT. 
Here and in the following, the subscript will always 
denote equilibrium quantities. 

The usual technique for simulating such a distribu- 
tion is to introduce a dynamics in configuration space. 
We choose a suitable set of transition rates, W[a — > a'\, 
which specify how a configuration a evolves into a new 
one, a', in unit time. For simplicity, we only consider 
transitions in which a and a 1 differ by a single nearest- 
neighbor particle-hole exchange. Now, the probability 
distribution, P(o~,t), becomes time-dependent and satis- 
fies a master equation (written, for simplicity, in contin- 
uous time): 



d t P(a,t) 



{W[a' -> cr]P(er', t) - W[cr -> a']P(a, t)} 

(2) 

Its stationary solution, P(cr) = limt—xx, P(cr, t), controls 
all time- independent properties. It is unique, under fairly 
generic conditions on the rates. To ensure that the de- 
sired equilibrium distribution, Po(<r), is reproduced, one 
chooses VF's which satisfy the detailed balance condition: 



W[a'^a] P (a) W 

Of course, this just ensures that every {...} bracket in Eq. 
(0) vanishes. An important quantity which enters here is 
the energy difference of two configurations, a' and a: 



A = H[a']-H[a] 



(4) 



A simple way of satisfying detailed balance is to impose 
rates which are functions of this difference alone: W[o~ — » 
a'] = w(/3Aq) where the function w must satisfy 



w(x) — w(—x)e 



(5) 



by virtue of Eq. Q but is otherwise arbitrary. All 
three rate functions to be considered below - Metropolis, 
Glauber, and heat bath - are constructed in this way, 
but differ in some important details. 

An obvious way of driving a system into a non- 
equilibrium steady state is to impose rates that violate 
detailed balance. A prototype model that has attracted 
much interest due to its remarkable properties is the 
driven Ising lattice gas (or KLS model) |3, La]- It dif- 
fers from the standard Ising model through the presence 
of an external force £, aligned with a particular lattice 
axis (the y-direction). When a particle attempts to jump 
to an empty nearest-neighbor site, it is not only affected 
by the local energetics, incorporated in Eq. |JTJ, but also 
by the drive: similar to an electric field, £ favors (sup- 
presses) particle hops along (against) the selected direc- 
tion, leaving transverse exchanges unaffected. A straight- 
forward extension of Eq. Q is to include the work done 
by the field, i.e., to define a local 'energy' difference of 
the form 



A = H[a']-H[a] - e£ 



(6) 
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Here, e = for two configurations differing only by a 
transverse jump, and e = +1 (— 1) if the particle hops 
along (against) the field in the move. We can now choose 
rates of the form JSJ with x — (3A. However, it is essen- 
tial to note that the combination of uniform drive and 
periodic boundary conditions precludes the existence of 
a global Hamiltonian for the driven system. A unique 
steady state, P(o-), establishes itself but cannot be ex- 
pressed in terms of a Boltzmann factor. To maximize 
the non-equilibrium effects, we choose infinite £ for our 
simulations, i.e., a particle will never jump against the 
field. 



Three different rate functions. 
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parallel to the drive. We also 



In this subsection, we introduce the three choices of 
transition rates - Metropolis, Glauber, and heat bath - 
which will be compared in the following. For the first 
two choices, the relevant quantity is the local 'energy' 
difference between the final (a') and the initial (a) con- 
figuration. For the third choice, the rate is independent 
of the initial configuration; instead, the selection crite- 
rion involves the local energy difference of the two pos- 
sible final configurations, a and a' . For the equilibrium 
Ising model, energy differences are easily computed from 
Eq. and each rate satisfies the detailed balance con- 
dition; for the driven model, we invoke Eq. ©, and de- 
tailed balance is violated. Random numbers are selected 
uniformly from the interval (0, 1). 

Metropolis dynamics. For this choice of rates, we ran- 
domly select a nearest-neighbor pair i, j of sites with 
different occupanices, o\ 7^ <7j. We denote the original 
configuration as a and let a 1 be the configuration with 
<7i, Oj switched (i.e., a- = ay <rj = o~i). The transition 
from a to a' is controlled by the Metropolis rate function, 
u>Met(x) = min(l,e~ 2: ). To be specific, we first compute 
x — j3A. If x < 0, the attempt (exchange) is accepted; 
if, however, x > 0, we draw a random number z and per- 
form the exchange only if z < e~ x . Clearly, energetically 
favorable moves are always performed while only a frac- 
tion of costly ones is accepted. As temperature increases, 
this fraction approaches 1 in a monotonic fashion. 

Glauber dynamics. Similar to Metropolis dynamics, 
the implementation of this algorithm involves, first, se- 
lecting two nearest-neighbor sites with different occu- 
pancies. Again, a' refers to the configuration with 
switched occupancies. Exchanges are then controlled by 
the Glauber rate function, wgi(x) = 1/(1 + e x ). Again, 
we compute x = (3A and draw a random number, z. If 
z < 1/(1 + e x ), we accept the exchange; otherwise, it 
is rejected. While energetically favorable moves are not 
necessarily accepted, they are always more probable than 
unfavorable ones. 

Heat bath dynamics. As pointed out above, the inter- 
pretation of a and a' is different here: These refer to the 
two possible final configurations of the central particle- 
hole pair. Showing only its local neighborhood in the 



(9) 



i=l 



i=4 



In equilibrium, the heat bath algorithm is of course 
isotropic: For both types of bonds, we select configura- 
tion a if a random number z satisfies z < 1/(1 + e~ 2/3h ); 
otherwise, we choose a' . For the driven case, this rule is 
only applied to bonds transverse to the drive; for paral- 
lel bonds, at infinite £, we choose configuration a with 
probability 1. 



(a) 



(b) 



(c) 



(d) 



FIG. 1: A central pair and a particular configuration of its 
six nearest neighbors. Occupied sites are indicated by solid 
circles. 

To appreciate the commonalities and differences of the 
three algorithms, it is useful to consider a simple example 
with infinite drive. Fig. 1 shows a central pair and a 
particular configuration of its six nearest neighbors. If 
the pair is aligned with the field direction (Figs, la, b), 
all three dynamics generate the same outcome: each will 
select Fig. la as the final configuration with probability 
1 for any value of (3 J. 

This is not the case for bonds transverse to the field 
(Figs, lc, d). For the purposes of this argument, we 
choose (3 J = 0.1. We denote the configuration shown 
in Fig. lc (d) by a (a 1 ). The 'energy' difference A = 
12 J is easily computed from Eq. (0J or @. Given a 
random number z, the Metropolis algorithm will accept 
a transition from a to a' only if z < e~ 12l3J ~ 0.30, while 
the reverse transition (a 1 to a) is always accepted. For 
Glauber dynamics, the transition from a to a' is accepted 
if z < 1/(1 + e 12/3J ) ~ 0.23, while the reverse transition 
is accepted if z < 1/(1 + e _12/3J ) ~ 0.77. Finally, the 
heat bath algorithm will choose a as the final state if z 
< 1/(1 + e~ 12l3J ) ~ 0.77, and a' otherwise. 
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The notable differences are these: First, the Metropolis 
algorithm accepts unfavorable moves with higher proba- 
bility than either heat bath or Glauber: e~ x > 1/(1 + e x ). 
As a result, it is more likely to explore unphysical do- 
mains of configuration space. Yet, it also accepts fa- 
vorable moves with higher probability, and thus leads 
to a more active dynamics. Comparing heat bath and 
Glauber rates, we note that both subdivide the unit in- 
terval into the same subsections (0.23 vs 0.77). Hence, 
they generate statistically very similar trajectories in con- 
figuration space. Update by update, however, the trajec- 
tories can differ: if, e.g., the initial configuration is a and 
the random number turns out to be 0.1, the heat bath 
algorithm will choose a' as the final configuration, while 
the Glauber rule leads to an exchange since z < 0.23. 
Yet, we will see below that this subtle difference does 
not affect the data. 



C. Structure factors and two-point correlations. 

Below their critical temperatures, both the Ising lat- 
tice gas and the KLS model phase-segregate into regions 
of high and low density, by virtue of the conservation 
law on the number of particles. Typical low-temperature 
configurations, for both models, show a single strip of 
high-density phase and its low-density mirror image. For 
the Ising model, the strip orients itself such as to mini- 
mize the energetic cost of interfacial length. In contrast, 
the low-temperature strip of the KLS model is always 
aligned with the direction of the drive, and the mini- 
mization of interfacial length does not play a dominant 
role (cf. Fig. 2). A quantity which easily distinguishes 
disordered configurations from such inhomogeneous ones 
is the (equal-time) structure factor. Written in spin lan- 
guage, it is defined as 



E 



e ik Vj 



(10) 



Here, k is a wave vector, taking discrete values k = 
(2Trn x /M, 2Tm y /L) with n x = 0, 1,...,M— 1 and n y — 
0, 1, ...,L — 1. For simplicity, we write S(n x ,n y ) in the 
following, and use 5(1,0) and 5(0,1) to detect strips 
aligned with the y- or cc-axis, respectively. For a per- 
fectly ordered strip aligned with y, 5(1,0) ~ 0A1ML 
is maximized; in contrast, a disordered configuration re- 
sults in 5(1,0) = 0(1). Further, the structure factor is 
the Fourier transform of the two-point correlation func- 
tion, G(r), defined via 



G (r) = (a a r ) - (ct ) (oy) 



(11) 



We assume translational invariance (modulo the lattice 
size) and invoke the half- filling constraint, whence (a r ) = 
((To) = 0. The same constraint imposes the sum rule 
J2 r G (r) = 5(0) = 0. Hence, negative values of G (r) 
for certain values of the argument should not come as a 
surprise. 



D. Energy histograms. 

For both the equilibrium Ising model and its driven 
counterpart, it is straightforward to accumulate a (nor- 
malized) energy histogram H{E,0), with respect to the 
energy function defined in Eq. For the equilibrium 
Ising model, H (E, (3) is intimately related to the Boltz- 
mann distribution: if W (E) denotes the density of states 
and Z(f3) the canonical partition function, we have 



H (E,f3) =Z- 1 (0)W(E) e 



-0E 



(12) 



Clearly, the right hand side is the probability Po(E,f3) 
to find the system with energy E. The power of the 
histogram method [12j resides in the observation that, 
up to statistical errors, a single histogram measured at 
temperature 1//3 is sufficient to construct Po(E, (3') at all 
other temperatures 1/(3': 



Po(E,P') = 



H (E,f3)e 



-(f3'-0)E 



J2 E ,H (E>,{3)e 



-{I3'-P)E> 



(13) 



This allows us to compute the moments of the energy 
distribution as functions of temperature and extract a 
wealth of thermodynamic information. 

For the driven lattice gas, H (E, (3) is easily compiled in 
a simulation. However, Eq. I|12|l certainly does not hold 
for the non-equilibrium steady state. In particular, exact 
solutions of small systems [l3| demonstrate unambigu- 
ously that, at a given temperature, configurations with 
the same energy need not have the same probability. At 
best, we can write, using the Kronecker symbol, 

H (E, P) = J2 5 e,h\*]P{o) = exp[-F(E, (3)} (14) 



where F(E, (3) is an as yet unknown function of its vari- 
ables which will certainly depend on the chosen dynam- 
ics. 

In the following, we probe F(E,(3) by considering a 
very simple ratio: We measure two histograms at differ- 
ent inverse temperatures, f3i and f3 2 and construct 



R(E, E') 



H{E, (3i) H{E',(3 2 ) 



H(E',/h) H(E,p 2 ) 



(15) 



for a range of E, E' . In equilibrium, this ratio is just a 
simple exponential: R (E,E') = exp[— I3 2 ){E— E')], 
since all normalization factors cancel. For the driven sys- 
tem, little is known except 

R{E,E') = ex V [F(E,p 2 )-F(E,p 1 )-F(E',^)+F(E',p 1 ) 

(16) 

This form will be analyzed further below. 

To conclude this section, we establish a few conven- 
tions and summarize the technical details of the simu- 
lations. All temperatures in the following are quoted 
in units of J/fcs; an important reference point is the 
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Onsager temperature T = -2/m(\/2 - 1) ~ 2.269 [3 
which marks the critical point of the two-dimensional 
Ising model. The equilibrium lattice gas and the driven 
system differ only in one parameter: £ = vs £ = 1000, 
respectively. Such a large value for £ suppresses (almost) 
all moves against the drive, and is therefore effectively in- 
finite. When a quantity, e.g., the critical temperature for 
the driven system, has been measured in different dy- 
namics, we will use superscripts M (Metropolis), H (heat 
bath), and G (Glauber) to distinguish them, as in T C M , 
T**, and T C G . The data for structure factors and two- 



point correlations were obtained on 100 x 100 systems 
while the histogram simulations used a smaller system 
size, 40 x 40. In each case, 1 MCS corresponds to one 
update attempt per site on average. For the larger sys- 
tem, each run lasted 2 x 10 6 MCS. The first 10 6 MCS 
were discarded to ensure that the system had reached 
steady state, and data were taken every 100 MCS over 
the second half of the run. For better statistics, 20 in- 
dependent runs were performed and averaged. For the 
smaller size, 4 x 10 6 MCS were discarded, followed by 
12 x 10 6 measurements. 



III. RESULTS 




FIG. 2: Typical configurations on a 48 x 432 lattice for the equilibrium Ising model using heat bath dynamics at T\ — 2.00 (a) 
and T 2 = 2.80 (b), and for its driven cousin at three temperatures: T 3 = 2.90 (c, d), T 4 = 3.30 (e, f) and T 5 = 3.70 (g, h). The 
first (second) configuration at each temperature was obtained using Metropolis (heat bath) dynamics. In each simulation, the 
data were collected after discarding 2 x 10 7 MCS for the equilibrium system and 10 7 MCS for its driven counterpart. 



A. Typical configurations, two-point correlations 
and structure factors. 



We begin our discussion by showing a few typical con- 
figurations of the driven system on a 48 x 432 lattice. 
Figs. 2a and b are obtained for the equilibrium case, just 
below and above criticality. The preference for horizontal 
interfaces is clearly seen in Fig. 2a. The remaining config- 
urations (Fig. 2c-h) all show the driven system, for heat 
bath and Metropolis dynamics, at three different temper- 
atures. At the lowest temperature T\ — 2.90 (Figs. 2c, d), 
the driven system is ordered for both dynamics. In stark 



contrast to the equilibrium case, the interfaces between 
high- and low-density regions are parallel to £ and there- 
fore clearly not dominated by energetics. At a slightly 
higher temperature, T2 = 3.30 (Figs. 2e, f), we observe 
the first glaring discrepancy between the two dynamics: 
the configuration generated by the heat bath algorithm 
is still ordered while the Metropolis configuration is al- 
ready disordered! Eventually, at T 3 = 3.70, both algo- 
rithms generate disordered configurations. Clearly, the 
two algorithms lead to different critical temperatures, 
with T C M < T]r . A rough estimate based on our data 
[H results in T C H = 3.55 ± 0.05 and T C M = 3.15 ± 0.05. 
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FIG. 3: The pair correlation function for the equilibrium system and its driven counterpart on a 100 x 100 lattice. The top row 
shows the Ising lattice gas at T = 2.47 with Metropolis (a), heat bath (b), and Glauber (c) dynamics; the bottom row shows 
the driven system at T = 3.60 with Metropolis (d), heat bath (e), and Glauber (f) dynamics. 
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FIG. 4: Parallel and transverse two-point correlations for the equilibrium system at T = 2.47 (a, b) and for the driven case 
at T = 3.60 (c, d) on a 100 x 100 lattice. Each plot shows data for three dynamics: Metropolis (asterisks), heat bath (filled 
squares) and Glauber (open circles). In (a, b) all data collapse, while in (c, d) only heat bath and Glauber data overlap. 



More precise estimates are available for Metropolis 
rates only: T C M = 3. 19801(19). 

To probe this apparent discrepancy between Metropo- 
lis and heat bath rates further, and to explore the posi- 
tion of Glauber rates in this triad, we turn to a more 



detailed analysis. In Fig. 3, we show surface plots 
of G (r) for the Ising lattice gas (top row, Figs. 3a-c) 
and the driven system (bottom row, Figs. 3d-f). The 
three columns correspond to the three different dynam- 
ics: Metropolis (Figs. 3a, d), heat bath (Figs. 3b, e), 
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FIG. 5: Structure factor contour plots for the equilibrium system and its driven counterpart on a 100 x 100 lattice. The top 
row shows the Ising lattice gas at T = 2.47 with Metropolis (a), heat bath (b), and Glauber (c) dynamics; the bottom row 
shows the driven system at T — 3.60 with Metropolis (d), heat bath (e), and Glauber (f) dynamics. 
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FIG. 6: Parallel and transverse structure factors for the equilibrium system at T = 2.47 (a, b) and for the driven case at 
T = 3.60 (c, d), on a 100 x 100 lattice. Each plot shows data for three dynamics: Metropolis (asterisks), heat bath (filled 
squares) and Glauber (open circles). Within error bars (not shown), the data effectively collapse in (a, b), while only heat bath 
and Glauber data overlap in (c, d). Note the different scale in (d). 



and Glauber (Figs. 3c, f). Fig. 4 shows selected projec- 
tions of G(r), namely G(0,y) and G(x,Q), for equilib- 
rium (Fig. 4a, b) and with infinite drive (Fig. 4c, d). As 
dictated by detailed balance, the correlation functions 
for the equilibrium system are independent of dynamics: 



there are no discernable differences between Figs. 3a-c, 
and the data in Figs. 4a and b collapse within statistical 
error bars (less than 0.01 in absolute units). The chosen 
temperature, T = 2.47, is close enough to Ising criti- 
cality so that lattice anisotropies are irrelevant: G(r) is 
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isotropic, with circular contours centered on the origin. 
The small negative values observed at large distances are 
a consequence of the sum rule. 

This simple picture becomes considerably more com- 
plex when we turn to the driven system (lower row 
of Fig. 3 and Figs. 4c, d). The chosen temperature, 
T = 3.60, is very close to our estimate for the critical 
temperature of heat bath and Glauber rates, T C H ~ T G 
~ 3.55 and about 15% above T C M . We immediately 
note the strong anisotropy induced by the drive. Fur- 
ther, there are noticeable differences between Metropo- 
lis rates on one hand, and heat bath and Glauber dy- 
namics on the other. These are most easily observed in 
Figs. 4c and d. For Metropolis rates, G M (0,y) is pos- 
itive and decreases monotonically throughout (Fig. 4c), 
while G M (x, 0) drops rapidly below zero, displays a min- 
imum and then recovers and approaches zero from be- 
low. These features have been noted before @, and 
are directly related to the breaking of detailed balance 
0, 0, . The data for Glauber and heat bath dynam- 
ics, while practically indistinguishable from one another, 
differ visibly from the Metropolis ones. Considering cor- 
relations measured along the field direction first, we ob- 
serve G H (0, y) ~ G G (0, y) > G M (0, y) for all y. In other 
words, Metropolis rates generate weaker correlations, 
consistent with the lower T C M . Heat bath and Glauber 
rates produce roughly the same correlations; moreover, 
these show clear signatures of being very close to criti- 
cality, evidenced by the distinctly positive value at the 
largest y shown: G ff (0,40) ~ G G (0,40) ~ 0.07. Highly 
correlated domains in the driven system are needle- 
shaped, with the needle pointing along the field, and this 
small, yet nonzero value indicates that some of these do- 
mains are long enough to span half the system. These 
precursors of ordering become even more obvious when 
we turn to correlations transverse to the field: The sec- 
ondary maximum in Fig. 4d indicates a tendency towards 
forming thin stripes for heat bath and Glauber rates. 

The structure factors bear out this picture. Again, 
the independence from the rates, and the isotropy near 
criticality is clearly displayed by the contour plots for 
the equilibrium system, shown in Figs. 5a-c, and by the 
projections shown in Figs. 6a, b. In the driven case 
(Figs. 5d-f and 6c, d), the presence of strong anisotropy 
is apparent, and the well-known discontinuity singularity 
at the origin [Tg is observed easily: lim^^+o S(k x , 0) 7^ 
linifc^o S(0, k y ). While these broad features characterize 
all three dynamics, the absolute values of the structure 
factors differ slightly from one another: S M (k x ,k y ) is 
generally smaller than either S H (k x ,k y ) or S G (k x ,k y ). 
Moreover, the distance from criticality can be measured 
through the discontinuity ratio, 

_ lim fcx ^ 5'(fc a; ,0) 
lim ky ^o S(0,k y ) 

which diverges as T — > T c Our data result in 

S M ~ 7.5, while S H ~ S G '~ 34 > S M . Our find- 



ings confirm, once again, that the heat bath and Glauber 
data are effectively much closer to criticality than those 
for Metropolis rates. 



0.02 F 




A = E - E' 

FIG. 7: Normalized histograms for the equilibrium system us- 
ing heat bath (solid line) and Metropolis (dotted line) rates, 
at /3i = 1/2.269 (a) and p 2 = 1/2.369 (b). In (c), we show 
ln.Ro vs E — E': Data are shown as open (Metropolis) and 
filled (heat bath) squares; the solid line is the expected be- 
havior, -(Pi -p 2 )(E-E'). 



B. Histogram Ratio Analysis 

In the final section, we turn to a brief investigation 
of energy histograms. Since Glauber and heat bath rates 
produce essentially identical data, we restrict ourselves in 
the following to just heat bath and Metropolis rates. To 
set the scene, we first show two histograms for the equi- 
librium system, generated at, and slightly above, crit- 
icality: T t = 2.269 and T 2 = 2.369 (Figs. 7a and b, 
respectively). As expected, the data for the different dy- 
namics collapse very well, within statistical errors. Not 
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FIG. 8: Normalized histograms for the driven system using 
heat bath (solid line) and Metropolis (dotted line) rates, at 
Pi = 1/3.550 (a) and p 2 = 1/3.650 (b). In (c), we show ]nR x 
vs E — E' . Metropolis data are shown as open circles and are 
taken at pi = 1/3.200 and fa = 1/3.300; the heat bath data 
(filled circles) are taken at p[ = 1/3.550 and P' 2 = 1/3.650. 
Two theoretical lines are shown: —(fix — P'2){E — E') (dotted) 
and -{fi[ -P' 2 )(E- E') (solid). 



surprisingly, the peak position shifts to higher energies 
with increasing temperature, while the width is largest 
at criticality. In Fig. 7c, we plot the corresponding his- 
togram ratio, Eq. (JTSJ, and compare it to the predicted 
exponential form. The agreement is of course very good. 

With Fig. 8, we enter novel territory. In analogy to 
the equilibrium plots, Figs. 8a and b display the energy 
histograms of the driven system, at two temperatures, 
Ti = 3.550 and T 2 = 3.650, for heat bath and Metropolis 
dynamics. The chosen temperatures correspond to criti- 
cality and slightly above for heat bath rates; for Metropo- 
lis rates, both are well inside the disordered phase. In 
contrast to the equilibrium case, the histograms clearly 



depend on the choice of rates: the peak positions are con- 
siderably higher for Metropolis than for heat bath rates. 
At the same time, the width is largest for the system clos- 
est to criticality, i.e., T\ = 3.550 with heat bath rates. 

A comment is in order, concerning the judicious choice 
of the two temperatures which enter the histogram ratio. 
It applies to both the equilibrium and the driven case. As 
we can see from Figs. 7 and 8, each histogram displays 
a well-developed peak. Energies far away from the peak 
position occur rarely, so that histograms are plagued by 
large statistical errors in those regions. In order to yield a 
reliable ratio, the corresponding histograms should over- 
lap in their statistically meaningful domains. Hence, the 
two chosen temperatures must not lie too far apart. 

In Fig. 8c, we present the histogram ratio for the driven 
system. For each dynamics, two temperatures close to 
their respective critical temperatures were chosen: 3.200 
and 3.300 for Metropolis rates, and 3.550 and 3.650 for 
heat bath rates. Remarkably, we observe that the his- 
togram ratio for both is again a simple exponential, i.e., 
lnR ao (E, E') oc {E — E'), at least over the range shown. 
In stark contrast to the equilibrium case, there is no a 
priori reason here to expect such behavior. Instead, it 
indicates that F(E, /3) in Eq. I|16fl depends sufficiently 
smoothly on E as to allow an expansion in A = E — E': 



In-RooGE, E') 



-A 
-a A 



9F(E,Pi) dF(E,(3 2 ) 



dE 

0(A 2 ) 



dE 



0(A 2 ) 
(17) 



Hence, the slope of the data in Fig. 8c allows us to 
probe a, as a function of temperature and dynamics. It 
manifestly differs from the equilibrium form — /3a)< A 
more systematic study is required to extract, and inter- 
pret, its properties. 



IV. CONCLUSIONS 

We have simulated the equilibrium Ising lattice gas and 
its driven non-equilibrium counterpart, using three differ- 
ent dynamics: Metropolis, Glauber and heat bath. In the 
equilibrium case, all three rate functions satisfy detailed 
balance with respect to the Ising Hamiltonian; as a con- 
sequence, all stationary (time-independent) equilibrium 
quantities are expected to be independent of the choice of 
the dynamics. Apart from unavoidable statistical errors, 
our equilibrium data are of course perfectly consistent 
with this expectation. For the driven system, this is no 
longer the case: due to the drive, all three rate functions 
violate detailed balance, and the 'decoupling' of station- 
ary properties from the chosen dynamics no longer holds. 
Measuring two-point correlations and structure factors in 
the disordered phase, we observe distinct differences be- 
tween the three dynamics. On the one hand, Metropolis 
rates lead to a lower critical temperature, and hence gen- 
erally weaker correlations in the disordered phase, than 
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either Glauber or heat bath rates. On the other hand, 
the latter two generate practically indistinguishable data. 
These features can be understood in terms of a few basic 
properties of the three rates: At a given temperature, 
Metropolis rates tend to accept all moves with a some- 
what higher probability than the other two rate func- 
tions. In other words, a system which evolves under the 
Metropolis algorithm 'sees' an effectively higher temper- 
ature than if it were running under heat bath or Glauber. 
A similar observation was made recently for field-driven 
Ising or solid-on-solid interfaces, subject to Glauber and 
Metropolis dynamics: there, Metropolis rates appear to 
lead to rougher interfaces and higher propagation veloc- 
ities than Glauber rates In contrast, heat bath and 
Glauber rates partition the unit interval into the same 
subsections and accept /reject moves according to this 
partition. As a result, they generate statistically indis- 
tinguishable trajectories in configuration space, leading 
to essentially identical data. 

It is essential to note, however, that the broad char- 
acteristics, associated with the breaking of detailed bal- 
ance, are clearly observed in all three dynamics: all struc- 
ture factors show the typical discontinuity singularity at 
the origin which, in turn, translates into power law de- 
cays of the two-point correlation functions. To summa- 
rize, universal features, associated with global symme- 
tries, remain independent of local changes of dynamic 
rules, both near and far from equilibrium. 

In a second part of this paper, we discuss the energy 



histograms H(E, (3) associated with our two models, gen- 
erated by heat bath and Metropolis dynamics. For the 
equilibrium system, the independence of the choice of 
dynamics is borne out again, while differences emerge in 
the driven case. A simple ratio, R(E,E'), constructed 
from two histograms measured at different temperatures 
[3\ and [3%, allows us to probe their functional form 
for a specified dynamics. In equilibrium, the canoni- 
cal distribution prescribes a simple exponential depen- 
dence, lni? = —(fti — th){E - E'). Remarkably, its 
non-equilibrium counterpart In is also exponential in 
(E — E'). This behavior indicates a smooth dependence 
of F(E, P) = — \n.H(E, (3) on E, allowing us to linearize 
In Rao in [E — E'). The slope of the resulting straight 
line depends on the dynamics and probes the derivative 
(dF/dE). Further, and more detailed, studies of this 
type may reveal some of the hidden "thermodynamics" 
of this remarkably complex non-equilibrium steady state. 
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